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ABSTRACT 

Context. The solar rotation profile is conical rather than cylindrical as one could expect from classical rotating fluid dynamics 
(e.g. Taylor-Proudman theorem). Thermal coupling to the tachocline, baroclinic effects and latitudinal transport of heat have 
been advocated to explain this peculiar state of rotation. 

Aims. To test the validity of thermal wind balance in the solar convection zone using helioseismic inversions for both the angular 
velocity and fiuctuations in entropy and temperature. 

Methods. Entropy and temperature fluctuations obtained from 3-D hydrodynamical numerical simulations of the solar convection 
zone are compared with solar profiles obtained from helioseismic inversions. 

Results. The temperature and entropy fiuctuations in 3-D numerical simulations have smaller amplitude in the bulk of the solar 
convection zone than those found from seismic inversions. Seismic inversion find variations of temperature from about 1 K at 
the surface up to 100 K at the base of the convection zone while in 3-D simulations they are of order 10 K throughout the 
convection zone up to 0.96 Rq. In 3-D simulations, baroclinic effects are found to be important to tilt the isocontours of f2 away 
from a cylindrical profile in most of the convection zone helped by Reynolds and viscous stresses at some locations. By contrast 
the baroclinic effect inverted by helioseismology are much larger than what is required to yield the observed angular velocity 
profile. 

Conclusions. The solar convection does not appear to be in strict thermal wind balance, Reynolds stresses must play a dominant 
role in setting not only the equatorial acceleration but also the observed conical angular velocity profile. 
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1. Introduction 

Helioseismic data from the Global Oscillation Network 
Group (GONG) and the Michelson Doppler Imager (MDI) 
have been used to infer the rotation profile in the solar 
interior (e.g., Thompson et al. 1996; Schou et al. 1998). 
The inversion results show that isocontours of the differ- 
ential rotation fl(r,9) are conical at mid-latitude, rather 
than cylindrical as was expected from early numerical 
simulations (e.g., Glatzmaier & Oilman 1982; Oilman & 
Miller 1986). More recent theoretical work (Durney 1999; 
Kitchatinov & Rudiger 1995; Brun & Toomre 2002 (here- 
after BT02); Rempel 2005; Miesch et al. 2006 (hereafter 
MBT06); Brun & Rempel 2008; Balbus et al. 2009) in- 
dicate that in order to break the Taylor-Proudman con- 
straint of cylindrical fl, the Sun must either have a system- 
atic latitudinal heat transfer in its convection zone or ther- 
mal forcing from the tachocline or most likely both. This 
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is due to the so-called thermal wind balance (Pedlosky 
1987), i.e., the existence in the solar convection zone of lat- 
itudinal entropy (or temperature) variation due to baro- 
clinic effect can result in a rotation state that breaks the 
Taylor-Proudman constraint. Such latitudinal variations 
of the thermal properties at the solar surface have been 
looked for observationally by several groups since the late 
60's (e.g., Dicke & Goldenberg 1967; Altroch & Canfield 
1972; Koutchmy et al. 1977; Kuhn et al. 1985, 1998, Rast 
et al. 2008; to cite only a few). This is a difficult task since 
one has to correct for limb darkening effect, photospheric 
magnetic activity, instrument bias and many other sub- 
tle effects to extract a relatively weak signal (see Rast et 
al. 2008). In most cases a temperature contrast of a few 
degree K is found from equator to pole at the surface, 
the pole being warmer. In some observations a minimum 
at mid latitude with a warm equator and hotter polar re- 
gions is also found. The warm polar regions and cool equa- 
torial region pattern is also found in 3-D simulation of the 
solar convection zone with temperature variation slightly 
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larger (i.e., of order 10 K; BT02, MBT06). At the surface a 
banded structure of the temperature field (warm-cool-hot) 
is also found in 3-D simulation of global scale convection. 
While very useful and instructive, most observations are 
confined to the solar surface and lack the information on 
the deep thermal structure of the solar convection zone 
which is key to characterise the dynamics of the deep so- 
lar convection zone. One way to remedy that limitation is 
to rely on helioseismic inversions that allows us to probe 
deeper in the Sun and to use 3-D global simulations of the 
solar convection zone to guide our physical understanding. 

Indeed, helioseismic inversions can give us the rota- 
tion rate, as well as the sound speed and density in the 
solar interior as a function of radius and latitude. Inside 
the convection zone the chemical composition is uniform 
and if we know the equation of state it is possible to 
determine other thermodynamic quantities like the tem- 
perature and entropy from the sound speed and density. 
Although, there may be some uncertainty in the equation 
of state, the OPAL equation of state (Rogers et al. 1996; 
Rogers & Nayafonov 2002) is quite close to the equation 
of state of solar material (e.g., Basu & Antia 1995; Basu 
& Christensen-Dalsgaard 1997). Thus, in this work we use 
the OPAL equation of state to calculate the perturbations 
in entropy and temperature and assess how well a strict 
thermal wind balance is established in the solar convec- 
tive envelope. To achieve this goal we make use of 2-D 
inversions of fi, S, T, using the GONG and MDI data for 
the full solar cycle 23 and analyse our findings using 3- 
D simulations obtained with the ASH (anelastic spherical 
harmonic) code (BT02; MBT06; Miesch et al. 2008) sup- 
ported by theoretical considerations on the thermal wind 
balance and vorticity equations. 

The rest of the paper is organised as follows: in Sect. 2 
we describe the data and technique used in this work while 
the results for the temperature and entropy inversions are 
described in § 3 along with those of 3-D simulations. In 
§ 4 we discuss at length the thermal wind balance and 
its generalisation and interpret our seismic inversion with 
3-D simulation of global scale convection. Finally, in § 5 
we put our results in perspective and conclude. 



2. The helioseismic data and inversion technique 

We use data from GONG (Hill et al. 1996) and SOI/MDI 
(Schou 1999). Eac;li data set consists of mean frequencies 
of different {n,l) multiplets, and the corresponding split- 
ting coefficients. We use 130 temporally overlapping data 
sets from GONG, each covering a period of 108 days, start- 
ing from 1995 May 7 and ending on 2008 May 9, with 
a spacing of 36 days between consecutive data sets. The 
MDI data consist of 61 non-overlapping data sets, each 
covering a period of 72 days, starting from 1996 May 1 
and ending on 2008 September 30. These data cover the 
solar cycle 23. For most of the work we use the temporal 
average over the available data to reduce the errors in in- 
version results. For this purpose we repeat the inversion 



process for all data sets and then take an average of all 
sets to get temporally averaged inversion results. 

We use a 2D Regularised Least Squares (RLS) in- 
version technique in the manner adopted by Antia et 
al. (1998) to infer the angular velocity in the solar in- 
terior from each of the available data sets. Similarly, we 
use a 2D RLS inversion technique as described by Antia et 
al. (2001) to infer the sound speed and density in the so- 
lar interior. In practice, we calculate the differences ScP /c^ 
and Sp/p with respect to a reference solar model. We use 
the solar model from Brun et al. (2002) with tachocline 
mixing as the reference model. In this work, we are only 
interested in the latitudinal variation in solar structure 
inside the convection zone. Thus the fiuctuation in sound 
speed can be converted to either temperature or entropy 
using the relation: 
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Here S is the specific entropy, T is temperature, P is 
pressure and Fi is the adiabatic index. The required par- 
tial derivatives are calculated using the OPAL equation 
of state. The derivatives of Fi are small in most of the 
convection zone, except for the ionisation zones of hydro- 
gen and helium, but for completeness we have included 
these derivatives in all our calculations. Helioseismic in- 
versions for rotation and asphericity are only sensitive to 
the North-South symmetric components and hence the in- 
verted profiles always show this symmetry. Hence, in this 
work we show the inversion results in only one hemisphere. 
Actual profiles may have some asymmetry about the equa- 
tor. 

3. Thermal perturbations in the solar convection 

zone 

Convection is a macroscopic transport of heat and energy. 
It is directly associated to correlation between the veloc- 
ity field and temperature fluctuations (Brun & Rempel 
2008). Being able to infer the temperature and entropy 
perturbations in the solar convection zone is thus key to 
understanding its turbulent dynamics. 



3.1. The Inverted profiles 

The aspherical part of temperature and entropy pertur- 
bations determined from temporally averaged GONG and 
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Fig. 1. The aspherical component of temperature fluctuation, 6T obtained from the temporally averaged GONG (left 
panel) and MDI (middle panel) data. The right panel shows the cuts at constant latitude of 5T obtained from MDI 
data along with la error estimates shown by dotted lines. All curves appear to merge at r = i?0, because 5T is of 
order of 1 K in that region. 






Fig. 2. The aspherical component of entropy fluctuation, 5S obtained from the temporally averaged GONG (left 
panel) and MDI (middle panel) data. The right panel shows the cuts at constant latitude of 5S obtained from MDI 
data along with la error estimates shown by dotted lines. 



MDI data are shown in Figures [T] and [21 The maximum 
temperature fluctuation near the bottom of the convec- 
tion zone is found to be about 100 K. These fluctuations 
increase with depth initially, because of steep increase in 
the temperature with depth which can induce an artifi- 
cially large value for 5T. The errors in 5T also increase 
with depth and the results may not be significant near 
the base of the convection zone. If we consider the rela- 
tive fluctuation 5T /T , then the maximum would be much 
closer to the surface and the value is of order of 10"'' or 
less. Similarly, if the entropy fluctuation are divided by its 
typical value of order Cp, then it too would be of the same 
order. Both these relative perturbations are of the same 
order as 5c^/c^. A detail look at Figures [T] and [5] reveal 
that the fluctuations are negative (relatively cold with re- 
spect to the spherically symmetric mean) at low latitude 
and warm at mid latitudes. In the bulk of the solar con- 
vection zone there is very little radial variation except 
near the surface. In the GONG data a cool polar region is 
also apparent but its significance is questionable given the 
relatively poor resolution of inversion techniques at high 
latitude. This feature is not clearly seen in the MDI data. 
While this latitudinal variation imprints through the sur- 



face for the entropy with little change in amplitude it is not 
the case for the temperature. At the surface the seismic in- 
version of the axisymmetric temperature fluctuations are 
very small in agreement with previous photospheric study 
(Rast et al. 2008). It may be noted that near the surface, 
above the lower turning point of the modes, the inversions 
may not be reliable. Around r = O.95i?0, where the inver- 
sions should be reliable, the temperature variations are of 
the order of 10 K. 

3.2. The profiles realised in 3-D models of large scale 
convection 

Recent efforts to develop high resolution global simula- 
tions of the solar convection zone in order to identify 
the physical processes at the origin of heat, energy and 
angular momentum transport have been quite successful 
at reproducing the seismically inverted differential rota- 
tion profile (BT02; MBT06). We display in Figured a 
typical solution of the solar convection zone and differ- 
ential rotation obtained with the ASH code (case AB3 of 
MBT06). Shown using a meridional cut is the longitudinal 
and temporal average of the angular velocity f2 along with 
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Fig. 3. Left panel: Angular velocity (in nHz) achieved in model AB3 of Miesch et al. 2006. Middle panel: Associated 
entropy S' (erg/g/K) and temperature T' (K) fluctuations with respect to the spherically averaged background. All 
quantities have been averaged over longitude and time (10 solar periods). Note the conical profile of the angular 
velocity at mid latitude and the latitudinal variation of the thermal variables possessing hot poles. Near the surface 
the temperature is banded with warm equator, cool mid latitudes and hot poles. Contrary to the helioseismic inversion 
no symmetry with respect to the equator is assumed and the quantities exhibit a small North-South asymmetry. 



the temperature and entropy fluctuations with respect to 
a spherically symmetric background. We first note that 
the differential rotation in the model is solar-like, with a 
fast equator and slow pole, and iso-contours of Q constant 
along radial lines at mid latitude (i.e, the rotation profile 
is conical rather than cylindrical). Its amplitude is also of 
the right order of magnitude. By contrast it is important 
to note that the temperature T' and entropy S' fluctu- 
ation^ are smaller by a factor of about 10 with respect 
to the seismic inversion, with temperature variations of 
about 10 K from equator to pole up to r = O.96i?0. A 
detailed analysis of the redistribution of heat and angular 
momentum in the 3-D models reveal that the Reynolds 
stresses and the latitudinal enthalpy flux are key players 
in establishing the profile of angular velocity and the vari- 
ation of temperature and entropy with latitude (Brun & 
Rempel 2008). Reynolds stresses transport angular mo- 
mentum from the polar region down to the equator being 
opposed by meridional circulation and viscous effect. The 
heat is transported poleward by the turbulent enthalpy 
flux (e.g. pCp{vgT'), with () denoting an azimuthal aver- 
age, p the mean background density and v'g the fluctuat- 

^ for the sake of clarity we make the distinction between the 
seimic inversion of the temperature and entropy perturbations 
denoted with a 5 symbol, and the one computed in the models 
denoted by a prime 



ing latitudinal component of the velocity field with respect 
to the axisymmetric mean, (see for more details Brun & 
Palacios 2009)) yielding cool equator and hot poles in most 
of the domain. It is opposed by thermal diffusion that tries 
to make the entropy and temperature field homogeneous. 
A careful study of the profile of temperature and entropy 
fluctuations reveals that the entropy is monotonic with 
respect to latitude while near the surface the temperature 
is banded (warm-cool-hot). Further the entropy profile is 
conical as is the angular velocity at mid latitude whereas 
the temperature profile is more cylindrical. In this strati- 
fied (anelastic) simulations the difference between the two 
thermal quantities is due to density (or pressure) fluctu- 
ations that cannot be neglected. This confirms that en- 
tropy is the key quantity to consider when studying the 
angular velocity profile of the Sun as clearly stated in the 
thermal wind equations detailed in §4.1. Mean field 2-D 
models also find axisymmetric temperature variations of 
order few Kelvin at the surface and in the bulk of the 
convection zone (Kitchatinov & Riidiger 1995, Kiiker & 
Riidiger 2005). Current global 3-D numerical simulations 
of the solar convection zone do not model the very surface, 
but stop at around 0.96 to 0.98 Rq, and as a consequence 
can not be used yet to model the near surface shear layer 
(see however the studies of Derosa et al. (2002) using a 



Brun, Antia & Chitre: Is the solar convection zone in strict thermal wind balance 



5 



modified ASH code or of Robinson & Clian (2001), using 
a spfierical wedge model). 

4. Quality of Thermal Wind Balance Achieved in 
the Sun and 3-D Models 

4.1. Theoretical Considerations 

In rotating convection, both radial and latitudinal heat 
transport occurs, establishing latitudinal gradients in tem- 
perature and entropy within the convective zone as illus- 
trated in Figure[31 A direct consequence of the existence of 
such gradients is that the surfaces of pressure and density 
fluctuations will not coincide anymore, thereby yielding 
baroclinic effects. We can turn to the vorticity equations 
( |Pedlosky(1987)[ |Zahn(1992)[ ) to analyse the role of the 
turbulence and baroclinic effects in setting the large scale 
flows shown in Figure [S] The thermal wind balance equa- 
tion can be derived from the vorticity equation as dis- 
cussed in detail by BT02 and MBT06. The equation for 
the vorticity in the purely hydrodynamical case can be 
derived under the anelastic approximation by taking the 
curl of the momentum equation (see also Fearn 1998 and 
Brun 2005 for its MHD generalisation and the notion of 
magnetic wind) : 



(wa • V)v - (v ■ V)a;„ - tVaiV ■ v) 
irVpxVP' -Vx i^er] - Vx(iv-X>), 



(4) 



'dt 



with U3a = Vxv -I- 2r2o the absolute vorticity, a; = 
Vxv the vorticity in the rotating frame, and V the vis- 
cous tensor given by: 
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(5) 



where is the strain rate tensor, and v is an effective 
kinematic viscosity. 

This vorticity equation helps in understanding the rel- 
ative importance of the different processes acting in the 
meridional planes. In the stationary case (^ = 0), and 
assuming an azimuthal average (such that ^ vanishes) 
the azimuthal component of Eq. ([4]) reads: 
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In the above equation we have identified several terms: 

• Stretching describes the stretching / tilting of the ab- 
solute vorticity due to velocity gradients; 

• Advection describes the advection of vorticity by the 
flow; 

• Compressibility describes the stretching of vorticity 
due to the flow compressibility; 

• ^ ^ is the baroclinic term, characteristic of non 
aligned density and pressure gradients; 

• -j— 4^ '^i^ ^ is part of the baroclinic term but arises 
from departure to adiabatic stratification; 

• Viscous accounts for the diffusion of vorticity due to 
viscous effects. 

We wish to stress that for the nonlinear stretching 
and advection terms (equivalent to Reynolds stresses 
in Navier-Stokes equation) their azimuthal average still 
yields partial derivatives in 0, since quadratic terms such 



^ r sin 9 dtp 



are non zero. 



Under the assumption that the convection zone is adi- 
abatic, the Rossby number Ro = lo /2VLq is small, and that 
compressibility, Reynolds and viscous stresses can be ne- 
glected, equation ([6]) simplifies to give: 



8z 



9 d{S') 
2V,orCp 89 



(8) 



Viscous stresses 



This is the thermal wind equation. It simply states that 
baroclinic effect can break Taylor-Proudman constraint of 
cylindrical differential rotation since otherwise 8vrf,/8z = 
for barotropic flows (Zahn 1992). This is due to the fact 
that the baroclinic terms drive meridional flows that under 
the influence of Coriolis force yield longitudinal flows that 
lead to a non cylindrical state of rotation. We now turn to 
our numerical simulation to evaluate the role played by all 
the terms of the vorticity equation identified above and to 
discuss the quality of the thermal wind balance achieved. 

4.2. Results from 3-D Models 

Figure[4]displays for case AB3 the left-hand side of Eq. ([6]) , 
along with the dominant terms of the right-hand side and 
their sum. We clearly see that the sum of the dominant 
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Fig. 4. Meridional cut of the terms discussed in equation 
(6) averaged over azimuth and 10 solar periods. Shown are 
in turn: d{v^)/dz, the stretching and advection of vortic- 
ity, the baroclinic effects, the viscous stresses and the sum 
of the RHS terms (we have divided all the RHS terms by 
l/2r2o)- Each panel share the same colour table with red 
denoting positive value. The min/max values used to scale 
the plots arc [—10^^, 10^^] expressed in . 



RHS term is in very close agreement with the LHS. We 
have chosen to form the temporal average over 10 solar pe- 
riods because it corresponds to about 10 convective over- 
turning times and lead to a very close balance between the 
LHS and the RHS of equation ([6|) . Shorter averages do not 
lead to such a good balance whereas longer averages do not 
change the quality of the balance obtained significantly 
nor the patterns of the various terms. Our more detailed 
decomposition of the vorticity equation is allowing us to 
identify which term is contributing and where. First the 
baroclinic term is found to be dominant in most of the 
bulk of the convection zone as was found by BT02 and 
MBT06. Advection terms are found to contribute both in 
the bulk and near the surface. They do not possess con- 
trary to the baroclinic term a systematic dominant con- 
tribution in each hemisphere. Their contribution leads to 
change in key places, yielding a more structured profiles 
of the RHS than the baroclinic term would have yielded if 
considered alone. Since the Rossby number realised in the 
simulation is less than one, we expect the advection and 
stretching term to be small on average in the simulation 
and indeed their maximum amplitude is not as large as the 
baroclinic term. However as stressed above this is not the 
case at all scales nor at all locations and they do contribute 
in key places, leading to the very good balance shown 
in Figure [4] between the LHS and RHS of equation ^ . 



Finally, in our models a viscous shear layer is dominat- 
ing the balance at the surface where the isocontours of VL 
possess the strongest latitudinal shear. Durney (1989) and 
Kitchatinov & Ruediger (1999) have also stressed that a 
strict thermal wind balance cannot be realised everywhere 
in the convection zone and that viscous stresses may play 
a role near the boundaries as observed in Figured We can 
thus conclude that equation ([8]) is only partly satisfied in 
our 3-D hydrodynamical simulations of the solar convec- 
tive envelope. Clearly baroclinic effects play a central role 
but these are far from being dominant everywhere and 
considering only equation ([5]) instead of the full balance 
expressed in equation ([6]), would be misleading. We now 
turn to seismic inversion to see if the thermal wind balance 
is strictly realised in the Sun or if other contribution must 
be invoked to explain the solar peculiar rotation profile. 

4.3. Inverted Solar Thermal Wind Balance 

The entropy perturbations obtained in section 13.11 can be 
differentiated to calculate the RHS of equation ([5]) 
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The result is shown in Figure \5[ We clearly see that the 
baroclinic term is non monotonic with respect to latitude, 
with large positive values near the poles and in a small re- 
gion at the equator whereas it is negative in mid latitudes. 
At the surface a surface thermal boundary layer is visible 
that yields strong radial gradients at high latitudes. 

As we have done with the 3-D model, the baroclinic 
term should be compared with 
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dr 
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which is shown in Figure [6l This quantity has much less 
structure in the bulk of the convection zone. Except for a 
slightly negative structure at mid depth and latitude, most 
of the action occurs in the surface shear layer where strong 
negative values arc found due to strong radial gradient of 
rotation rate in the near surface shear layer. This near 
surface layer is not present in the simulations and hence 
can not be compared with the results of 3-D simulations. 
It is clear that contrary to what we have seen with the 3-D 
model in the previous section, the two quantities do not 
agree with each other even slightly. In fact, these two terms 
differ by more than an order of magnitude. While the term 
involving in equations 8 or 10 is of the same order in 
both the 3-D simulations and the seismic inversion (^ 
10^^ s^^), this is not the case for the baroclinic terms due 
to the very large entropy and temperature variations in 
the inverted profiles. Although the profile in simulations 
qualitatively reproduces the features seen in solar profile, 
the detailed latitudinal variations in the two do not match 
precisely. 
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Fig. 5. The aspherical component of the latitudinal derivative of the entropy fluctuation, {g /2norCp){dSS/ dO) obtained 
from the temporally averaged GONG (left panel) and MDI (right panel) data. The values are in s~^. 





-2x10-' -1.5x10-' 



-IxlO-a -5xl0-' 



-2x10-' -1.5x10-' 



■IxlO-a -5x10-" 



Fig. 6. The derivative of rotation velocity, sm 9 {r cos 980,/ dr — sm9dQ,/dO), obtained from the temporally averaged 
GONG (left panel) and MDI (right panel) data. The values are in s^^ . 



5. Discussion of Results 

What can be the source of the disagreement between the 
inverted baroclinic contribution and the z derivative of 
the angular velocity (i.e equations 8, or 9 and 10)? The 
first and easiest solution is that the inversion of the ther- 
mal quantities lack the necessary accuracy and given the 
increase by two orders of magnitude of the background 
temperature and density with depth, we end up with vari- 
ations that are too large. The source of discrepancy will 
then be due to an overestimation of 5T and 5S. It is 
not easy to decide if these inverted thermal fluctuations 
are too large or if the simulations (both 2-D and 3-D) 



underestimates the fluctuations realised in the Sun, be- 
cause for instance of their limited Reynolds number. We 
must thus also consider the possibility that these large 
thermal perturbations are genuine. If this is indeed the 
case we need to see how we could resolve the discrepancy 
between the seismically inverted LHS and RHS of equa- 
tion 8. As stated in section 14.11 to obtain a strict ther- 
mal wind balance as expressed in equation 8, one need 
to make a certain number of assumptions: adiabaticity, 
weak Rossby number, negligible compressibility, viscous 
and Reynolds stresses, stationarity. Further by consider- 
ing only the hydrodynamic contributions we have omitted 
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those associated with Maxwell stresses that are certainly 
present in the magnetic Sun. Let's however assume that 
the Maxwell stresses are not the source of the large ob- 
served discrepancy. We are confident that this is the case 
because we have formed temporal averages over a maxi- 
mum and minimum period of activity and the differences 
between the two periods are about 10 times smaller that 
what it would required if all the sources of discrepancy 
were coming from the Maxwell stresses alone. We never- 
theless intend to make a more systematic study of the 
departure of the strict thermal wind balance linked to 
magnetic effects (i.e. via the so called magnetic wind) by 
analysing the solar cycle 23 in details and by comparing 
with dynamo simulations of the solar convection (Brun 
et al. 2004). We must thus question the validity of the 
other hypothesis made in deriving equation 8. Clearly it 
is justified given the very low microscopic value of the so- 
lar kinematic viscosity to consider that the viscous terms 
do not contribute much. This is clearly not the case in 
the 3-D models where near the surface they are major 
contributors to the overall balance (see Figure [H middle 
panel of the bottom row), but this is due to our large 
effective viscosity. Assuming adiabaticity is certainly rea- 
sonable in most of the convection zone but clearly not 
near the surface. Since we are mostly interested in under- 
standing the bulk dynamics of the solar convection zone, 
this term is indeed very small. The choice of low Rossby 
number that allows us to neglect uj over 2Qo the planetary 
vorticity is certainly not justified at all scales of the tur- 
bulent velocity spectra, in particular for those scales much 
smaller than the Rossby radius of deformation (Pedlosky 
1987). In the Sun the large range of convection scales cer- 
tainly undergo different dynamics depending on how sen- 
sitive they are to the Coriolis force. The subtle angular 
momentum and heat redistribution realised in the Sun is 
in part captured in our 3-D models. We can thus anal- 
yse if the Reynolds stresses associated with the turbulent 
motion indeed play a central role. As discussed in detail 
in Brun & Toomre (2002) and in §4 we know that it is 
indeed the case in our numerical simulations (see Figure 
m middle and right panel of the top row) even though 
our simulation do not possess a Reynolds number and a 
degree of turbulence as high as that in the Sun. We can 
thus expect, given the very large Reynolds number of the 
solar convection zone, that Reynolds stresses must play 
a central role in the Sun in shaping the differential rota- 
tion profile and that they somehow in part compensate 
the baroclinic contribution to yield the observed profile of 
angular velocity. This is a very interesting results, since 
it indicates that the differential rotation is of dynamical 
origin for its equatorial acceleration (as revealed by study- 
ing angular momentum transport in our simulations as in 
BT02 or Miesch et al. 2008) but also most certainly for its 
shape with Reynolds stresses helping or opposing in some 
regions the baroclinic effects to break Taylor-Proudman 
constraint. Of course this conclusion only holds if the in- 
verted large thermal fluctuations are real. 
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